package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.TP3;
import static scu.maqiang.numeric.Constants.TW3;
import static scu.maqiang.numeric.MVO.average;
import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES2T3PE extends FES {

	public FES2T3PE(Mesh mesh) {
		super(mesh, 3);
		eleStiff = FES2T3PE::elementPEStiff;
		eleMass = FES2T3PE::elementPEMass;
		eleSource = FES2T3PE::elementPESource;
	}

	/**
	 *
	 * @param xy
	 * @param mat
	 * @param tp
	 * @param K
	 */
	public static void elementPEStiff(int[] ele, double[][] xy, double[][] mat, BVPType tp, double[][] K) {
		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
		double[][] grad = new double[3][2];
		double area = CG.gradient2DT3(xy, grad);
		double[][] B = null;
		double r0 = average(x);
		MVO.fill(K, 0.0);
		switch (tp) {
//		case PLANESTRAIN:
//		case PLANESTRESS:
		case CONSTITUTE_MATRIX_COMMON:
			B = new double[5][9];
			B[0][0] = B[2][1] = B[3][2] = grad[0][0];
			B[0][3] = B[2][4] = B[3][5] = grad[1][0];
			B[0][6] = B[2][7] = B[3][8] = grad[2][0];
			B[1][1] = B[2][0] = B[4][2] = grad[0][1];
			B[1][4] = B[2][3] = B[4][5] = grad[1][1];
			B[1][7] = B[2][6] = B[4][8] = grad[2][1];
			mul(matmul(transpose(B), matmul(mat, B)), area, K);
			break;
		case CONSTITUTE_MATRIX_AXIS:
			B = new double[6][9];
			area = area * r0;
			B[0][0] = B[3][1] = B[4][2] = grad[0][0];
			B[0][3] = B[3][4] = B[4][5] = grad[1][0];
			B[0][6] = B[3][7] = B[4][8] = grad[2][0];
			B[2][0] = B[2][3] = B[2][6] = 1.0 / (3 * r0);
			B[1][1] = B[3][0] = B[5][2] = grad[0][1];
			B[1][4] = B[3][3] = B[5][5] = grad[1][1];
			B[1][7] = B[3][6] = B[5][8] = grad[2][1];
			mul(matmul(transpose(B), matmul(mat, B)), area, K);
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public static void elementPEMass(int[] ele, double[][] xy, double[][] rho, BVPType tp, double[][] K) {
		double[] p = rho[0];
		double area = CG.area2DT3(xy);
		double[] x = new double[] { xy[0][0], xy[1][0], xy[2][0] };
		if (tp == BVPType.CONSTITUTE_MATRIX_AXIS || tp == BVPType.AXIS_LUMP_MASS) {
			MVO.mul(p, x, p);
		}
		double N11 = (3 * p[0] + p[1] + p[2]) / 30 * area;
		double N22 = (p[0] + 3 * p[1] + p[2]) / 30 * area;
		double N33 = (p[0] + p[1] + 3 * p[2]) / 30 * area;
		double N12 = (2 * p[0] + 2 * p[1] + p[2]) / 60 * area;
		double N13 = (2 * p[0] + p[1] + 2 * p[2]) / 60 * area;
		double N23 = (p[0] + 2 * p[1] + 2 * p[2]) / 60 * area;
		MVO.fill(K, 0.0);
		switch (tp) {
		case CONSIST_MASS:
		case AXIS_CONSIST_MASS:
			K[0][0] = K[1][1] = N11;
			K[0][3] = K[1][4] = K[3][0] = K[4][1] = N12;
			K[0][6] = K[1][7] = K[6][0] = K[7][1] = N13;
			K[3][3] = K[4][4] = N22;
			K[3][6] = K[4][7] = K[6][3] = K[7][4] = N23;
			K[6][6] = K[7][7] = N33;
			break;
		case LUMP_MASS:
		case AXIS_LUMP_MASS:
			K[0][0] = K[1][1] = (2 * p[0] + p[1] + p[2]) / 12 * area;
			K[3][3] = K[4][4] = (p[0] + 2 * p[1] + p[2]) / 12 * area;
			K[6][6] = K[7][7] = (p[0] + p[1] + 2 * p[2]) / 12 * area;
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
		// System.out.println(MVO.toString(K));
		// System.exit(0);
	}
	
	public void assemblePESource(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS) {
		assembleGlobalVector(coefFunc, param, FES2T3PE::elementPESource, tp, RHS);
	}

	public void assemblePESource(double[] constCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(constCoef, FES2T3PE::elementPESource, tp, RHS);
	}

	public void assemblePESource(double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(vecCoef, FES2T3PE::elementPESource, tp, RHS);
	}
		
	public static void elementPESource(int[] ele, double[][] xy, double[][] vf, BVPType tp, double[] eleF) {
		double area = CG.area2DT3(xy);
		double[] fu = vf[0];
		double[] fv = vf[1];
		double[] rf = vf[2];
		double[] x = MVO.col(xy, 0);
		double sumFu = MVO.sum(fu);
		double sumFv = MVO.sum(fv);
		double sumRf = MVO.sum(rf);
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case CONSTITUTE_MATRIX_COMMON:
			for (int i = 0; i < 3; i++) {
				eleF[3 * i] = (fu[i] + sumFu) * area / 12.0;
				eleF[3 * i + 1] = (fv[i] + sumFv) * area / 12.0;
				eleF[3 * i + 2] = (rf[i] + sumRf) * area / 12.0;
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
		case AXISYMMETRIC:
			for (int j = 0; j < 4; j++) {
				for (int i = 0; i < 3; i++) {
					eleF[3 * i]     += area * TW3[j] * MVO.dot_product(TP3[j], fu) * MVO.dot_product(TP3[j], x) * TP3[j][i];
					eleF[3 * i + 1] += area * TW3[j] * MVO.dot_product(TP3[j], fv) * MVO.dot_product(TP3[j], x) * TP3[j][i];
					eleF[3 * i + 2] += area * TW3[j] * MVO.dot_product(TP3[j], rf) * MVO.dot_product(TP3[j], x) * TP3[j][i];
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}

	public void assembleSurfaceForce(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T3PE::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public void assembleSurfaceForce(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES2T3PE::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public void assembleSurfaceForce(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T3PE::elemenSurfaceForce, tp, RHS, bdLabel);
	}

	public static void elemenSurfaceForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double length = CG.length2DL2(xy);
		double[] gu = f[0];
		double[] gv = f[1];
		double sumGU = MVO.sum(gu);
		double sumGV = MVO.sum(gv);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case CONSTITUTE_MATRIX_COMMON:
			length /= 6.0;
			for (int i = 0; i < 2; i++) {
				eleF[3 * i] = (sumGU + gu[i]) * length;
				eleF[3 * i + 1] = (sumGV + gv[i]) * length;
			}
			break;
		case AXISYMMETRIC:
		case CONSTITUTE_MATRIX_AXIS:
			length /= 12.0;
			double[] x = MVO.col(xy, 0);
			double sumX = MVO.sum(x);
			for (int i = 0; i < 2; i++) {
				eleF[3 * i] = (sumGU * sumX + 2 * gu[i] * x[i]) * length;
				eleF[3 * i + 1] = (sumGV * sumX + 2 * gv[i] * x[i]) * length;
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
	
	public void assembleSurfaceElectricDisp(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS,
			int[] bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T3PE::elemenElectricDisp, tp, RHS, bdLabel);
	}

	public void assembleSurfaceElectricDisp(double[] constCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(constCoef, FES2T3PE::elemenElectricDisp, tp, RHS, bdLabel);
	}
	
	public void assembleSurfaceElectricDisp(double[][] vecCoef, BVPType tp, double[] RHS, int[] bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T3PE::elemenElectricDisp, tp, RHS, bdLabel);
	}

	public static void elemenElectricDisp(int[] ele, double[][] coord, double[][] g, BVPType tp, double[] eleF) {
		double len = CG.length2DL2(coord);
		double sumG = MVO.sum(g[0]);
		double sumX = coord[0][0] + coord[1][0];
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 2; i++) {
				eleF[3 * i + 2] = (sumG + g[0][i]) * len / 6.0;
			}
			break;
		case AXISYMMETRIC:
			for (int i = 0; i < 2; i++) {
				eleF[3 * i + 2] = (sumG * sumX + 2 * g[0][i] * coord[i][0]) * len / 12.0;
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
	}
}
